# R function for making a stacked coefficient plot from different regression estimates
# by Piero Stanig, based in part on Yu-Sung Su's coefplot code
# =====================================================================
# November 22, 2006
# This version: May 2021

stackplot<- function(fit,
                     longnames=NULL, center.labels=FALSE,return.position=FALSE,
                     x.label="", y.label="", main.label="", sub.label="",
                      varnames=TRUE, estimators=FALSE, one.sd=FALSE,fifty=TRUE,
				lwd.sd=1,max.n=15, offset=.1, smallish=.02,
                     box=FALSE, symmetric=TRUE, cex.var=1.2, cex.pts=1.5,  suppress.intercept=TRUE, thin.box=FALSE,
                     restrict=NULL, suppress=NULL, ind.labels=NULL, pch.type=21, smallmargin=FALSE
                     ){

 


		if (!(class(fit)=="matrix"|class(fit)=="data.frame"
		     )){

			stop(message = paste(class(fit), "class not supported!\n"))
		}
	




if (class(fit)[1]=="matrix"|class(fit)[1]=="data.frame"){

	coeffs<-as.matrix(fit)
	
	n.models<-ncol(coeffs)/2
	
	estimators<-FALSE

		
	if (is.null(longnames)){

		if (!is.null(row.names(fit))){

			longnames<-row.names(fit)
			
		}else longnames<-" "
	
	}

	n.vars<-nrow(coeffs)
	
	max.length<-n.vars

	if (max.length>40) {
	
		stop(message = "Can display up to 40 coefficients per model in a plot!\n")

	}

}

if (max.length>40) {

	stop(message = "Can display up to 40 coefficients per model in a plot!\n")

}


####Manage the variable names

if (is.null(longnames)){

	longnames<-predictor.names[1:max.length,]

	longnames<-t(na.exclude(t(longnames)))[,1]
}

####Stack the coefficients and the sd's

coef.mat<-array(NA, c(n.models*n.vars,2))

suppress.mat<-array(NA, c(n.models*n.vars,1))

	for (i in 1:n.models){
		
		for (j in 1:n.vars){

			coef.mat[((j-1)*n.models)+i,  1] <-coeffs[j,(i*2)-1]
      
      suppress.mat[((j-1)*n.models)+i,  1] <-is.element(i,suppress)

      
			coef.mat [i+(j-1)*n.models, 2]<-    coeffs[j,i*2]
      
      
		}
	}

coef<-coef.mat[,1]


sd<-coef.mat[,2]

omit.from.plot <- as.numeric(suppress.mat)
	#calculate the HPDs
	
missing.c<-is.na(coef)

missing.s<-is.na(sd)

coef[missing.c==TRUE]<-0

sd[missing.s==TRUE]<-0



###now stack them

n.x <- length(coef)

####GROUP THE idx SO THAT THE COEFFICIENTS 
####ARE CLOSE TO EACH OTHER

idx <- rep(seq(1 , n.x/n.models),each=n.models)

smallsteps<-rep(.4*seq(0,1,length.out=n.models), n.x/n.models) 

idx<-idx+smallsteps

min.x <- round(min(coef-2*sd), 1)-0.1           
                    
max.x <- round(max(coef+2*sd), 1)+0.1                                 

x.scale <- c(min.x-offset, max.x+offset)                                            

y.scale <- range(idx)                   

if(offset!=.1){symmetric=FALSE}
if (symmetric) x.scale<-c(-max(abs(min.x), abs(max.x)), max(abs(min.x), abs(max.x)) )



subtitle<-""

if (estimators){

	esti.names<-array(NA, c(1,n.models))

	for (k in 1:n.models){

	
	}
	subtitle<-esti.names[1,1]
	
		if (n.models>2){
			
			for (k in 2:(n.models-1)){

				subtitle<-paste(subtitle, esti.names[1,k], sep=", ")

			}
		 }
	
	if (n.models>1) subtitle<-paste(subtitle, esti.names[1,n.models], sep=" and ")
	
	subtitle<-paste(subtitle, "estimates")               
	
}

if (subtitle=="") subtitle<-sub.label

notshown <- as.numeric(missing.c)+as.numeric(omit.from.plot)
notshown[notshown==2]<-1

if(smallmargin==FALSE){
par(mar=c(4+1*as.numeric(estimators==FALSE), 
2, 2+1*as.numeric(main.label!=""), 
.2+.6*max(nchar(longnames)))   )
}else{ 
  
  par(mar=c(3.5+1*as.numeric(estimators==FALSE), 
            2, 2+1*as.numeric(main.label!=""), 
            .1+.22*max(nchar(longnames)))   )  }
               
plot(coef, idx, axes=FALSE, type="n",                                     
       
xlim=x.scale, ylim=y.scale+c(-.02, .02), 
        
xlab=x.label, ylab=y.label,                   
        
main=main.label, sub=subtitle)                                                        


axis(1)    
             
#axis(3)

abline(v=0, lty=2, col="grey40",lwd=1-.4*thin.box)                                                 

if(n.models==1|length(pch.type)==1){
  
points(coef, idx, pch=pch.type, cex=cex.pts*(1-notshown))

  } else{
    
    temp.i =1 
    
for( i.dots  in 1:length(idx)){
  
  points(coef[i.dots], idx[i.dots], pch=pch.type[temp.i], cex=cex.pts*(1-notshown[i.dots]))
  
  temp.i=temp.i+1 
  
  if(temp.i> length(pch.type)){temp.i=1}
  
}
  }  
if (box==TRUE){
	
	for (k in 2:(n.vars)){

	abline(h=k-.2, lty=3, lwd=1-.4*thin.box)
	}

}

  
  if(fifty==TRUE){crit=1}else{crit=1.645}
	if (one.sd==FALSE) segments ((1-notshown)*(coef+2*sd), 
	idx, (1-notshown)*(coef-2*sd), idx, col="Grey80", lwd=lwd.sd)


	segments ((1-notshown)*(coef+crit*sd), idx, 
            (1-notshown)*(coef-crit*sd), idx, lwd=lwd.sd)     
		
		  




if(center.labels==TRUE){adj.here=.5}else{adj.here=1}
if (varnames==TRUE) {axis(4, .2*as.numeric(n.models>1)+(1:n.x), 
		longnames[1:n.x], las=2, tck=TRUE, pos=1.1*max(coef+2*sd), lty=0, adj=adj.here, cex.axis=cex.var)
  
  out.now <- 1:n.x
}
  
		
if(!is.null(ind.labels)){
  
text(coef, idx-smallish, ind.labels, cex=.8)  
  
}		
par(mar=c(5, 4, 4, 2) + 0.1)


if(return.position==TRUE){
  
  return(out.now)}

}
